// sort matrix rows by TSP algorithm
// use pearson correlation distance
#include <iostream>
#include <fstream>
#include <boost/ptr_container/ptr_vector.hpp>
#include <boost/integer_traits.hpp>
#include <boost/property_map/property_map.hpp>
#include <boost/graph/adjacency_matrix.hpp>
#include <boost/graph/metric_tsp_approx.hpp>
#include <vector>
#include <numeric>
#include <boost/tuple/tuple.hpp>
#include <string>
#include "util.hh"

typedef std::vector<double> vec_t;
typedef boost::ptr_vector<vec_t> mat_t;

bool read_matrix_file(const char* mat_file, mat_t& mat)
{
  size_t n = num_rows(mat_file);
  size_t p = num_columns(mat_file);

  if (n < 2)       // must have at least two rows
    return false;

  mat.clear();
  double x;
  std::ifstream ifs(mat_file, std::ios::in);

  for (size_t j = 0; j < n * p && (ifs >> x); ++j)
  {
    size_t r = j / p;
    size_t c = j % p;
    if (c == 0)
      mat.push_back(new vec_t(p, 0.));
    mat[r][c] = x;
  }
  ifs.close();

  return true;
}

typedef boost::adjacency_matrix<boost::undirectedS, boost::no_property,
    boost::property<boost::edge_weight_t, double,
        boost::property<boost::edge_index_t, int> > > Graph;

template<class DistMeasure>
void construct_graph(const mat_t& dist_mat, Graph& g)
{
  size_t n = dist_mat.size();

  boost::property_map<Graph, boost::edge_weight_t>::type weight = boost::get(
      boost::edge_weight, g);

  Graph::edge_descriptor e;
  bool ok;
  DistMeasure dist;
  for (size_t i = 0; i < n; ++i)
  {
    for (size_t j = (i + 1); j < n; ++j)
    {
      double d = dist(dist_mat[i], dist_mat[j]);

      boost::tie(e, ok) = boost::add_edge(i, j, g);
      if (ok)
      {
        boost::put(weight, e, d);
        // std::cerr << i << ", " << j << " = " << d << std::endl;
      }
    }
  }
}

struct square_op_t
{
  double operator ()(double a, double x)
  {
    return a + x * x;
  }
};

struct euclidean_dist_t
{
  double operator ()(const vec_t& x, const vec_t& y)
  {
    double ret = 0.;
    for (size_t i = 0; i < x.size(); ++i)
    {
      double d = x.at(i) - y.at(i);
      ret += d * d;
    }
    return std::sqrt(ret);
  }
};

struct pearson_dist_t
{
  //      sum_i x_i * y_i - n * xbar * ybar
  // r = ------------------------------------------------------------------
  //      sqrt{ sum_i x_i^2 - n * xbar^2 } sqrt{ sum_i y_i^2 - n * ybar^2 }
  double operator ()(const vec_t& x, const vec_t& y)
  {
    assert_msg(x.size() == y.size(), "x and y must have same length");

    double n = x.size();
    double xsum = std::accumulate(x.begin(), x.end(), 0.);
    double ysum = std::accumulate(y.begin(), y.end(), 0.);
    double x2sum = std::accumulate(x.begin(), x.end(), 0., square_op_t());
    double y2sum = std::accumulate(y.begin(), y.end(), 0., square_op_t());

    double xy_sum = 0.;
    for (size_t i = 0; i < x.size(); ++i)
      xy_sum += x.at(i) * y.at(i);

    double r = (xy_sum - xsum * ysum / n) / std::sqrt(x2sum - xsum * xsum / n)
        / std::sqrt(y2sum - ysum * ysum / n);

    return 1. - r;
  }

};

int main(const int argc, const char* argv[])
{
  if (argc < 2)
  {
    std::cerr << argv[0] << " matrix file [metric]" << std::endl;
    std::cerr << "metric = pearson or euclidean" << std::endl;
    std::cerr << "default = pearson" << std::endl;
    return 1;
  }

  const char* mat_file = argv[1];
  mat_t dist_mat;

  assert_msg(read_matrix_file(mat_file, dist_mat), "cannot read matrix");

  Graph g(dist_mat.size());

  if( argc > 2 )
  {
    std::cerr << "metric = " << argv[2] << std::endl;
    if( argv[2] == "euclidean" )
      construct_graph<euclidean_dist_t>(dist_mat, g);
    else
      construct_graph<pearson_dist_t>(dist_mat, g);
  }
  else
  {
    construct_graph<pearson_dist_t>(dist_mat, g);
  }

  typedef Graph::vertex_descriptor Vertex;

  std::vector<Vertex> order;
  boost::metric_tsp_approx_tour(g, std::back_inserter(order));

  size_t n = dist_mat.size();
  assert_msg(order.size() == (n + 1), "must have a tour");

  for (size_t i = 0; i < n; ++i)
    std::cout << order[i] << std::endl;

  return 0;
}
